function [sigmaVec, covMat, snr, snrVec, condNum, condNumVec, Params]...
	 = computeerrbounds_network(Params, detNameVec, noisePSDVec, lambdaVec)
%
% COMPUTEERRBOUNDS_NETWORK - compute the RMS errors in estimating the parameters
% of a coalescing binary using a netowrk of detectors, employing 3.5PN SPA waveforms 
% 
% usage: [sigmaVec, covMat, snr, snrVec, condNum, condNumVec...
%           ] = computeerrorbounds(Params, detNameVec, noisePSDVec, lambdaVec)
%
% P. Ajith, 8 Aug 2009
% 
% $Id: computeerrbounds_network.m 68 2010-01-20 21:58:45Z anand.sengupta $

 setconstants

 % compute the PN errors, if appropriate 
 flso = (1./sqrt(6))^3 / (pi*Params.totalMass*MSOLAR_TIME);

 if flso > Params.fLower

    if length(detNameVec) ~= length(noisePSDVec)
	error('### length(detNameVec) is not equal to length(noisePSDVec)');
    end

    % compute the fisher matrix for each detector
    for iDet=1:length(detNameVec)

	% load the detector tensor
	Params.Det = [];
	Params.Det = LoadDetectorInfo(detNameVec{iDet}, Params.Det, Params.gmstRad);

	% Calculate the A and B (direction dependent parameters of the source) and P,Q
	[Params.A, Params.B] = CalculateAtBt(Params.Det, Params.Alpha, Params.Delta);

	% Calculate the P and Q parameters. Aeff = P + iQ. See Notes.
	if Params.isAngAvged
	   
	    Params.P = sqrt(35*(Params.A^2 + Params.B^2)/32)...
			* (Params.totalMass*MSOLAR_TIME)^(5/6) * pi^(-2/3)...
			/ Params.r;

	    Params.Q = -Params.P;

	else
	
	    Params.P = sqrt(5*Params.eta/96)...
	        	* (Params.totalMass*MSOLAR_TIME)^(5/6) * pi^(-2/3)...
		 	* (Params.a1*Params.B + Params.a3*Params.A);

	    Params.Q = -sqrt(5*Params.eta/96)...
		 	* (Params.totalMass*MSOLAR_TIME)^(5/6) * pi^(-2/3)...
		 	* (Params.a4*Params.B + Params.a2*Params.A);
	end

	% compute the fisher matrix
	[fishMatDet, snrVec(iDet), AeffWeightdDet]...
		 = computefishermatrix_ifo( Params, noisePSDVec{iDet}, lambdaVec);

	% condition number of the fisher matrix
	condNumVec(iDet) = cond(fishMatDet);

	if iDet == 1
		fishMat = fishMatDet;
		AeffWeightd = AeffWeightdDet; 
	else
		% sum the fisher matrices of the individual detectors
		fishMat = fishMat + fishMatDet;

		% sum the noise-weighted (fourier-domain) amplitude of the
		% signal at individual detectors -- in order to compute the 
		% network SNR
		AeffWeightd = AeffWeightd + AeffWeightdDet; 
    	end

    end

    rankMat = rank(fishMat);

    condNum = cond(fishMat);

    % if the condition number is sufficiently small, invert the fisher matrix
    % to compute the covariance matrix. If not, fill the covariance matrix 
    % with zeros. 
    if condNum < 1e30
	covMat = inv(fishMat);
    else
	warning('condition number > 1e30. setting covariance matrix to zero');
	covMat = zeros(size(fishMat));
    end

    % compute the RMS errors from the diagonal elements of the covariance matrix
    sigmaVec = sqrt(diag(covMat))';

    % compute the network SNR
    snr = 2*sqrt(sum(AeffWeightd.^2))

 else 
	warning('flso is greater than fLower');
	sigmaVec 	= [];
	covMat 		= [];
	snr 		= [];
	snrVec 		= []; 
	condNum 	= [];
	condNumVec 	= [];
 end

 return;


